function [longTermSeq, guideFunction] = GenerateLongTermSequence(seqLen,period_labels)

if nargin < 3
    numTargets = 7;
end

if nargin < 2
    period_labels = 7:15;
end

if nargin < 1
    seqLen = 128;
end


if mod(seqLen,2) ~= 0
    error('seqLen should be even');
end



% Generate guide function
guideFunction = zeros(1,seqLen);

for thisPeriod = period_labels
    guideFunction = guideFunction + sin((2*pi/thisPeriod)*(1:seqLen) + (rand * pi));
end

guideFunction = guideFunction - min(guideFunction);
guideFunction = guideFunction ./ max(guideFunction);


% Generate sequence
longTermSeq = zeros(1,seqLen);

for i=1:seqLen
    longTermSeq(i) = binornd(1,guideFunction(i));
end


% Adjust the sequence so that it has the same number of 1s and 0s
excessOfOnes = sum(longTermSeq) - (seqLen/2);
[~,sequenceOfIndicesToChange] = sort(abs(guideFunction-0.5));

while excessOfOnes ~= 0
    
    index_to_change = 1;
    % If there are more 1s than 0s, remove some 1s
    if excessOfOnes > 0
        % Find the preferred index to flip
        while longTermSeq(sequenceOfIndicesToChange(index_to_change)) == 0
            index_to_change = index_to_change + 1;
        end
        % Flip it!
        longTermSeq(sequenceOfIndicesToChange(index_to_change)) = 0;
        
    else % If there are more 0s than 1s, remove some 0s
        % Find the preferred index to flip
        while longTermSeq(sequenceOfIndicesToChange(index_to_change)) == 1
            index_to_change = index_to_change + 1;
        end
        % Flip it!
        longTermSeq(sequenceOfIndicesToChange(index_to_change)) = 1;
        
    end
    
    excessOfOnes = sum(longTermSeq) - (seqLen/2);
end


% Shift the sequence, so that is has values 1 and 2
longTermSeq = longTermSeq + 1;


%plot(longTermSeq);hold on; plot(guideFunction,'r');axis([1 128 -0.1 1.1]);

